subroutine calculate_characteristics_normal(points,vertex2D,nVer2D,nInCat2D,&
    params,numbersEl,time,MN)
    common /characteristics/ nPoints,nMn
    common /dimTask2D/ nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    common /commonInfo/ eps2d,eps3d,radiusRED,maxiter,typeTask,nt,nMaterials
    integer elemsCount2D,typeTask,timeLevel,count
    character*4 Ext
    real*8 eps2d,eps3d,radiusRED,time,vertex2D,params,points,strength_normal,x,y,r,z,MN,&
    xM,yM,zM,xN,yN,zN,dx,dy,dz,Ex
    dimension time(nt),vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),nInCat2D(elemsCount2D),&
    params(nMaterials2d),points(3,nPoints),numbersEl(4,nodesCount2D),strength_normal(3),&
    MN(6,nMn),Ex(nt)

    if(nPoints.ne.0) then
        do j=1,nPoints
            call searchElementForPoint(points(1,j),points(2,j),vertex2D,nVer2D,&
            numbersEl,nEl1,nEl2)
            write(Ext,'(i4.4)') j
            open(500,file='.\results\SP_E_NX.'//Ext,mode='write')
            write(500,"(/,A)") '      t (mc)         ExN (mB)'
            open(501,file='.\results\SP_E_NY.'//Ext,mode='write')
            write(501,"(/,A)") '      t (mc)         EyN (mB)'

            do timeLevel=10,nt
                r = dsqrt(points(1,j)**2+points(2,j)**2)
                if(r.eq.0.d0)r=1d-3
                z = 0.d0
                call calculateStrength(vertex2D,nVer2D,params,nInCat2D,nEl1,nEl2,&
                points(1,j),points(2,j),r,z,strength_normal,timeLevel)
                write(500,*)time(timeLevel)*1000.d0,strength_normal(1)*1000.d0
                write(501,*)time(timeLevel)*1000.d0,strength_normal(2)*1000.d0
            enddo

            close(500)
            close(501)
        enddo
    endif

    if(nMn.ne.0) then
        count = 20
        do j=1,nMn
            xM = MN(1,j)
            yM = MN(2,j)
            zM = MN(3,j)
            xN = MN(4,j)
            yN = MN(5,j)
            zN = MN(6,j)
            dx = (xN-xM)/count
            dy = (yN-yM)/count
            dz = (zN-zM)/count
            do timeLevel=10,nt
                Ex(timeLevel) = 0.d0
            enddo

            do k=1,count
                x = xM+dx/2.d0+(j-1)*dx
                y = yM
                call searchElementForPoint(x,y,vertex2D,nVer2D,numbersEl,nEl1,nEl2)
                do timeLevel=10,nt
                    r = dsqrt(x**2+y**2)
                    if(r.eq.0.d0)r=1d-3
                    z = 0.d0
                    call calculateStrength(vertex2D,nVer2D,params,nInCat2D,nEl1,nEl2,&
                    x,y,r,z,strength_normal,timeLevel)
                    Ex(timeLevel) = Ex(timeLevel) + strength_normal(1)*dx
                enddo
            enddo

            write(Ext,'(i4.4)') j
            open(500,file='.\results\SP_dV_NX.'//Ext,mode='write')
            write(500,"(/,A)") '      t (mc)         dV(mV)'
            do timeLevel=10,nt
                write(500,*)time(timeLevel)*1000.d0,Ex(timeLevel)*1000.d0
            enddo
            close(500)
        enddo
    endif
end
!*************************************************************************
subroutine searchElementForPoint(x,y,vertex,nVer,numbersEl,nEl1,nEl2)
    common /dimTask2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    real*8 x,y,vertex,r,z
    integer elemsCount2D
    dimension vertex(2,nodesCount2D),nVer(4,elemsCount2D),numbersEl(4,nodesCount2D)

    r = dsqrt(x**2+y**2)
    if(r.eq.0.d0) then
        r=1d-3
    endif
    z = 0.d0

    nEl1 = elemsCount2D
    do while(nEl1.ge.1.and.&
             (r.lt.vertex(1,nVer(1,nEl1)).or.&
              r.gt.vertex(1,nVer(2,nEl1)).or.&
              z.lt.vertex(2,nVer(1,nEl1)).or.&
              z.gt.vertex(2,nVer(3,nEl1))))
        nEl1 = nEl1-1
    enddo

    nEl2=numbersEl(3,nVer(1,nEl1))
end
!*************************************************************************
subroutine calculateStrength(vertex,nVer,params,nInCat2D,nEl1,nEl2,x,y,r,z,&
    strength_normal,timeLevel)
    common /dimTask2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,nodesTermCount2D,&
    m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    common /commonInfo/ eps2d,eps3d,radiusRED,maxiter,typeTask,nt,nMaterials
    real*8 x,y,vertex,r,z,eN,hz1,hz2,hz3,params,sigma1,sigma2,ValueInPoint2d,&
    strength_normal,Hfi1,Hfi2,Hfi3,Hfi4,Hfi5,Hfi6
    integer elemsCount2D,timeLevel
    character*4 Ext
    dimension vertex(2,nodesCount2D),nVer(4,elemsCount2D),eN(4),params(nMaterials2d),&
    nInCat2D(elemsCount2D),strength_normal(3)

    write(Ext,'(i4.4)') timeLevel-1
    open(502,file='.\calculations\Hfi\v2.'//Ext,access='direct',recl=8,mode='read')
    read(502,rec=nVer(1,nEl1))Hfi1
    read(502,rec=nVer(2,nEl1))Hfi2
    read(502,rec=nVer(3,nEl1))Hfi3
    read(502,rec=nVer(4,nEl1))Hfi4
    read(502,rec=nVer(1,nEl2))Hfi5
    read(502,rec=nVer(2,nEl2))Hfi6

    hz1                = vertex(2,nVer(1,nEl1))-vertex(2,nVer(1,nEl2))
    hz2                = vertex(2,nVer(3,nEl1))-vertex(2,nVer(1,nEl1))
    sigma1             = params(nInCat2D(nEl1))
    sigma2             = params(nInCat2D(nEl2))
    k1                 = hz1/(hz1+hz2)
    k2                 = hz2/(hz1+hz2)

    eN(1)              = -(Hfi3-Hfi1)/sigma2/hz2*k2-(Hfi1-Hfi5)/sigma1/hz1*k1
    eN(2)              = -(Hfi4-Hfi2)/sigma2/hz2*k2-(Hfi2-Hfi6)/sigma1/hz1*k1
    eN(3)              = -(Hfi3-Hfi1)/sigma2/hz2
    eN(4)              = -(Hfi4-Hfi2)/sigma2/hz2

    strength_normal(1) = ValueInPoint2d(eN,vertex,nVer,nEl1,r,z)*x/r
    strength_normal(2) = ValueInPoint2d(eN,vertex,nVer,nEl1,r,z)*y/r
    close(502)
end
!*************************************************************************
real*8 function ValueInPoint2d(v,vertex,nVer,num,r,z)
    common /dimTask2D/nodesCount2D,elemsCount2D,nodesBoundCount2D,&
    nodesTermCount2D,m2D,mTerm2d,nodesNotTermCount2D,nMaterials2d
    integer elemsCount2D
    real*8 v,vertex,psi1,psi2,psi3,psi4,r1,r2,z1,z2,r,z
    dimension v(4),vertex(2,nodesCount2D),nVer(4,elemsCount2D)
    r1 = vertex(1,nVer(1,num))
    r2 = vertex(1,nVer(2,num))
    z1 = vertex(2,nVer(1,num))
    z2 = vertex(2,nVer(3,num))
    ValueInPoint2d = psi1(r1,r2,z1,z2,r,z)*v(1)+&
                     psi2(r1,r2,z1,z2,r,z)*v(2)+&
                     psi3(r1,r2,z1,z2,r,z)*v(3)+&
                     psi4(r1,r2,z1,z2,r,z)*v(4)
end
!*************************************************************************
!*************************************************************************
!*************************************************************************